#
"""
 赵振华 2023-04-06
 用scipy.integrate 积分
 傅里叶级数展开
"""
#
import numpy as np
import matplotlib.pyplot as plt
import scipy as sci

def f(x):  # 注意修改函数类型记得修改函数的积分区间
    # return np.sin(x) * x + 1  # 目标函数
    # return np.cos(x)**3
    # return np.sin(x)**3
    return np.sin(x)


def g(x):  # 原函数的级数，展开到 number_k 阶
    y = [a0] * len(x)
    for j in range(1, number_k):
        y = y + a_k[j] * np.cos(2 * np.pi * j * x / periodic) + \
                b_k[j] * np.sin( 2 * np.pi * j * x / periodic)
    return y
 
low = 0# -np.pi  # 周期的下界
high = np.pi  # 周期的上界
periodic = high - low  # 周期
x = np.linspace(low, high, 50)  # 定义采样点的个数
number_k = 50  # 展开到 number_k 阶
# 求展开系数 a_k， b_k
a0 = sci.integrate.quad(f, low, high)[0] / periodic  # 展开系数
a_k = np.zeros(number_k+1)
b_k = np.zeros(number_k+1)
for j in range(0, number_k):
    if j == 0:
        a_k[0] = a0
        b_k[0] = 0 # b_0 = 0
    else:
        a_k[j] = 2 *sci.integrate.quad(lambda x: f(x) * np.cos(2 * np.pi * j * x / periodic), low, high  )[0] / periodic
        b_k[j] = 2 *sci.integrate.quad(lambda x: f(x) * np.sin(2 * np.pi * j * x / periodic), low, high  )[0] / periodic


plt.figure(1)
plt.plot(x, f(x), color="red",label='Original function Sin(x)')
plt.scatter(x, g(x), color="blue", linewidth=1,label='Fourier series')
plt.legend()
plt.savefig('fourier_series_plot1.png', dpi=300)
plt.show()

plt.figure(2)
k = np.linspace(0, number_k, number_k+1)
plt.scatter(k, a_k, color="red",label='a_k')
plt.scatter(k, b_k, color="blue",label='b_k')
plt.legend()
plt.savefig('fourier_series_plot2.png', dpi=300)
plt.show()
